A study of a local Monte Carlo technique for simulating systems of charged particles 
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We study some aspects of a Monte Carlo method invented by Maggs and Rossetto for simulating 
systems of charged particles. It has the feature that the discretized electric field is updated locally 
when charges move. Results of simulations of the two dimensional one-component plasma are 
presented. Highly accurate results can be obtained very efficiently using this lattice method over 
a large temperature range. The method differs from global methods in having additional degrees 
of freedom which leads to the question of how a faster method can result. We argue that efficient 
sampling depends on charge mobility and find that the mobility is close to maximum for a low rate 
of independent plaquette updates for intermediate temperatures. We present a simple model to 
account for this behavior. We also report on the role of uniform electric field sampling using this 
method. 

PACS numbers: 87.10.+e,87.15.Aa,41.20.Cv 
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I. INTRODUCTION 

As an aid to understanding the behavior of statistical 
systems, it is common practice to carry out computer 
simulations of models that are invented to capture the 
physical behavior of interest because, even for unelab- 
orate models, it is often extremely difficult to find ap- 
proximate and especially exact analytical results. Both 
Monte Carlo and molecular dynamics methods work by 
generating large numbers of particle or field configura- 
tions from which correlation functions can be obtained 
that converge to thermal averages by design. In Monte 
Carlo, each new configuration is produced by making a 
pseudorandom change to the last one in the chain, testing 
the ratio of Boltzmann weights (which, in the canonical 
ensemble, amounts to a computation of the energy differ- 
ence), and choosing either the new configuration or the 
old one. Molecular dynamics is a numerical evolution 
of a many-body system from the equations of motion of 
the component 'particles' with some constraint depend- 
ing on the ensemble - for example to maintain a constant 
temperature. The degree of approximation to thermal 
equilibrium in both cases is limited, in practice, by the 
efficiency with which the algorithm samples the regions of 
highest probability in the space of configurations, and the 
rate at which configurations can be produced. These two 
factors determine the range of system sizes that can be 
sampled effectively, though a compromise must often be 
made between obtaining large numbers of configurations 
based on the relaxation time of the combined model and 
algorithm and choosing large enough systems to lessen 
finite size effects that obscure the thermal equilibrium re- 
sults required and which are usually poorly understood. 

This work is a study of a Monte Carlo technique for 
simulating systems of thermalized classical charges. Most 
of the systems studied in soft condensed matter labo- 
ratories fall into this class. It includes simple liquids, 
polymers, colloids, liquid crystals, proteins and other bi- 
ologically active molecules. 

In many interesting cases, charged particle interactions 



can be screened by the charges themselves, giving an ef- 
fective interaction of the Yukawa form with a screening 
length that can often be tuned in experiments. But there 
are also many cases where one must consider the exact 
Coulomb interaction. Because the Coulomb interaction 
is long-ranged, one must presumably either compute all 
pairwise interactions with the charge that is moved at a 
given step (in Monte Carlo) or one must solve the Poisson 
equation at each step. Consider a Monte Carlo simula- 
tion with N charges. When a charge is moved, the en- 
ergy difference between configurations requires O(N) op- 
erations for all pairwise interactions. Each Monte Carlo 
sweep, which is defined to be N charge moves or Monte 
Carlo steps, requires 0(N 2 ) operations; this is how the 
simulation will slow down as the system size increases. 
With short range interactions, in contrast, we expect 
O(N) operations for each sweep. 

The order 0(N 2 ) scaling for Coulomb interactions is a 
real problem and a great deal of effort has been devoted 
to improving the scaling and the prefactor that accom- 
panies it. This study focuses on a relatively new method 
proposed by Maggs and Rossetto [l|, which achieves 
0(N) scaling and which has some promise as a com- 
petitive approach for simulating charged particles using 
either molecular dynamics or Monte Carlo. 



A. Local Electrostatics 

The gist of the method of Maggs and Rossetto is as fol- 
lows. Rather than solving the Poisson equation (which 
follows from Gauss' law supplemented with the condi- 
tion that the electric field circulation vanish) at each 
step, one imposes only Gauss' law. The electric field 
circulation decouples from the charges, so it averages out 
during the simulation and one obtains electrostatic av- 
erages. Gauss' law can be imposed locally in the sense 
that when a charge moves, this divergence condition can 
be maintained by changing the electric field in the neigh- 
borhood of the charge only. Transposing this result to 
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Monte Carlo, we see that each sweep requires only O(N) 
operations. 

Now we fill in details of the argument [l| . Consider a 
number of charged particles with charge density p(r) = 
J2i1i^( r ~~ r i)- Suppose the electric field E is affected 
by the charge distribution only to the extent that Gauss' 
law is satisfied 

V • E = p(r). 

In other words we do not impose the further condition 
required for electrostatics that the curl of the electric field 
be equal to zero. The charged particle system is closed so 
that the amount of charge does not change and the whole 
system is immersed in a heat bath at inverse temperature 
[3. The partition function in this canonical ensemble in 
two dimensions is 

Z((3) = jY[d 2 r t ypE5(V-E-p(r)) 

i 

xexp(^ J d 2 rE 2 (r)), (1) 

in which the factor resulting from the momentum inte- 
gration over each particle has been left out. The electric 
field can be split into an irrotational or longitudinal part 
E|| and a transverse part E^: E = Ey + E^ satisfying 

V • E_l = and V x E|| = 0. 

Gauss's law can now be written as V • E|| = p(r) and 
with the condition on the longitudinal field that its curl 
vanish, E|| is derivable from a scalar field 0(r) which 
obeys Poisson's equation 

V 2 0(r) = -p(r) 
E|, = -V#r). 

The solution to the Poisson equation in two dimensions 
is uniquely 

0(r)=-^- / d 2 rV(r')ln|r-r'|. 

27T J 

We return to the energy functional to see the effect of 
separating out the longitudinal and transverse contribu- 
tions 

J d 2 r E 2 (r) = J d 2 r [(V0) 2 + E^ - 2V0 • Ej_] . 

Partially integrating the third term, we get 

2 J d 2 r<j){\7 -EjJ = 0. 

The whole of the dependence of the energy on the charge 
positions is contained in the scalar potential so the par- 
tition function factorizes: 

z(p) = z E± (p) >< / n^ 2 ^ ex p (-§ / ^^(w) 2 ) 



where Z E± (0) , which comes from the Gaussian integra- 
tion over the transverse degrees of freedom, is indepen- 
dent of the charge configuration. So thermodynamic av- 
erages derived from the partition function |T]) are the 
same as those that would be found from the partition 
function with the electrostatic energy functional. 

Because Gauss' law can be imposed locally, we ought 
to be able to run Monte Carlo simulations of systems 
of charged particles with local field updates and hence 
local energy computations at each step rather than solv- 
ing Poisson's equation to find the new unique global field 
configuration for each charge configuration. The advan- 
tage of having local energy computations at each step 
is an improvement in the scaling of the algorithm over 
methods involving direct evaluation of the pairwise inter- 
action. Each particle move is a local operation involving 
0(1) operations, so the scaling of the time per Monte 
Carlo sweep is 0(N). The transverse field must also be 
sampled, and yet one finds that this does not contribute 
significantly to the scaling. This latter fact is quite coun- 
terintuitive and we shall devote some more attention to 
it in the following. 

It is worth pointing out that if we begin with Maxwell's 
electrodynamics, a similar argument to the one given 
above shows that electrostatic averages follow in this 
case. This is not remarkable given the decoupling of 
the divergence and circulation of the electric field in the 
pure electric field case described above, but we set out 
this result in an appendix. One could therefore devise a 
molecular dynamics algorithm with the appropriate con- 
straint for working at constant temperature that inte- 
grates Maxwell's equations (which are local equations) 
at each time-step and with the speed of light as a vari- 
able, chosen to be small enough so that the algorithm is 
truly local, and obtain the same results that would be 
obtained directly from electrostatics. However, Maggs 
and Rossetto have shown that one is free to find the sim- 
plest local dynamics possible for the charges that simul- 
taneously samples longitudinal and transverse degrees of 
freedom subject to Gauss' law, and use this dynamics to 
perform much more efficient Monte Carlo or molecular 
dynamics simulations than would be possible with full 
Maxwell electrodynamics. In this paper, we shall use a 
charge and field dynamics invented by Rottler and Maggs 
Q that omits magnetic fields altogether. 



B. Overview 

In the next section, we provide a very short introduc- 
tion to charge interactions in periodic media in the con- 
tinuum limit and for the case where the electric field is 
discretized onto a lattice. Hence, we describe the conven- 
tional 0{N 2 ) global Monte Carlo method we have used 
and give a practical account of the corrections that must 
be made to lattice simulations to better approximate 
the continuum pair- wise potential. Section IIIII describes 
the implementation of the local O(N) method for Monte 
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Carlo simulations with charges moving off-lattice with 
particular attention to unusual features of the dynamics; 
I also present a continuum version of the local update 
method which shows some already known features very 
transparently. Global and local Monte Carlo methods are 
applied in Section llVBI to the 2D OCP and we discuss the 
accuracy of the local methods and the effect of uniform 
field sampling. Then, we shall see that charge autocorre- 
lation times are much greater for the local method than 
for global methods. Section H V D I includes the autocorre- 
lation data into a comparison of the speeds of the local 
method and a global method for typical system sizes. Re- 
sults of an investigation into the effects of the transverse 
electric field on charge relaxation times are given in Sec- 
tion HVE] We also discuss the problem of how additional 
degrees of freedom result in O(N) scaling when the naive 
expectation is that transverse fields should be allowed to 
relax across the system at each step. The first appendix 
derives a useful result for corrections to 2D lattice-based 
simulations. The second appendix shows that the argu- 
ment of Maggs and Rossetto given above carries over to 
Maxwell electrodynamics. 



II. SIMULATION METHODS 

The Monte Carlo simulations that will be discussed be- 
low are based on the acceptance procedure of Metropolis 
and co-workers This involves the computation of a 
difference in energies of two field configurations at each 
step. The simulation methods we have used are distin- 
guished by the ways in which this difference is found. 

In this section, we discuss the interactions between 
charged particles in rectangular cells with periodic 
boundary conditions and the corrections that have to be 
made to the interaction when the electric field is put onto 
a lattice in order to better approximate the continuum 
potential energy. We describe the simulation methods we 
have used in turn, starting with the continuum case, then 
the problems of discretization, and leaving the special 
properties of the local algorithm until the next section. 
In the following, unless otherwise stated, we consider two 
dimensional electrostatics for notational simplicity and 
because this is most appropriate to the simulations that 
we have performed. 



A. Periodic systems 

The potential energy of a finite number of charges in 
a rectangular cell with periodic boundary conditions is a 
conditionally convergent series over all pairwise Coulomb 
energies for the charges in the infinite array of cells im- 
plied by the boundary conditions. Each cell must be 
charge neutral otherwise the series is divergent. And if 
the dipole moment in each cell is zero the sum is ab- 
solutely convergent. In general, the dipole moment in 
each cell is non-vanishing, so the problem is to devise a 



physically motivated summation of the series. A further 
problem is to identify the unique solution to the Pois- 
son equation with periodic boundary conditions. These 
have been solved by Leeuw and Perram [J] by obtain- 
ing the potential for a large spherical cluster of identi- 
cal cells embedded in a medium of dielectric constant e; 
their expression includes a term that depends on the net 
dipole moment of the elementary cell and the dielectric 
constant of the surrounding medium. The unique peri- 
odic solution to the Poisson equation is obtained in the 
limit as e — > oo and corresponds to the case of a large 
spherical cluster of identical rectangular cells of charges 
surrounded by a perfect conductor that perfectly cancels 
the dipole moment of the cluster of cells The peri- 
odic potential, called the Ewald potential, f (|Al[) given 
in the first appendix for the two dimensional case [J]) 
is a pairwise potential for charges in a rectangular box 
that includes the direct Coulomb potential and also the 
potential due to all the image charges for the pair consid- 
ered. The dipole moment term that is present for finite 
e takes the form 

U ^=2LMe + l) (f> a ) (2) 

for cell lengths L x and L y and with charge q a at position 
r a . This will be discussed further in connection with the 
local method. 

The Ewald potential splits the calculation of the po- 
tential into short-range and long-range parts. The lat- 
ter is performed in Fourier space. This forms the ba- 
sis for molecular dynamics simulations. However, the 
Ewald potential is not always the most efficient starting 
point for computations. I have used the form found by 
Lekner [6|, |7| and developed for rectangular periodic cells 
by Gr0nbech- Jensen If the charges are separated by 
a vector with Cartesian components x and y and the cell 
lengths are, once again, L x and L y , then the pairwise 
potential energy between unit charges is 




where A' is a constant that includes the interaction of the 
charges with the uniform background charge and with all 
their image charges and is given by 

*(£)-i£-5"--n('M-<^))- 

The advantage of this formula is that it is in terms of 
elementary functions; although it is an infinite product 
in its main part, the product converges very quickly - I 
take \k\ < 5 which agrees with the exact result to machine 
precision. 
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B. Lattice simulations 

The O(N) Monte Carlo method that is the subject of 
this paper requires that we discretize the electric field so 
that it lives on the links of a square mesh with charges on 
the sites. Overall charge neutrality must be maintained 
(if necessary by having a uniform background charge). 
We consider the case of a two-dimensional square lat- 
tice with lattice spacing a, site labels {R = (n x , n y )} 
and link directions fi = x,y such that charge <?(R) is 
assigned (in the way described above) to site R and elec- 
tric fields are defined on the links of the lattice with 
E(R, —fx) — —E(TH — fi, fj.) by convention. The lattice 
has N x — L x /a, N y = L y /a plaquettes in the x and y 
directions respectively. Sites n x = and n x = N x are 
identified. Periodicity in the y direction is implemented 
in the same way. 

The discrete form of Gauss' law is 

Y, a (£( R > m) + £(R, -m)) = g(R). (4) 

When the Poisson equation is discretized onto the lat- 
tice, its solution can be written as 

0(R) = £]G(R-R')«(R')- 

R' 

The charge at lattice site R is q(R) and G(R) = 
G(n x ,n y ) is the lattice Green's function. One can show, 
by calculating the discrete Fourier transform of the Pois- 
son equation with periodic boundary conditions, that 

G(n X: n y ) = Y exp ( 2 ^^r^) E exp ( ^^Q 1 - ) 



N x N y 4 _2cos(^)-2cos(^) 

in two dimensions. The electric field is related to the 
potential on the lattice by aE x (n Xl n y ) — 4>{n x + 1, n y ) — 
(j)(n Xl n y ) and similarly for the y component. 

The charged particles are not confined to the lattice 
sites in our simulations. In common with most molecular 
dynamics simulations of charged particles, charges are 
smoothly assigned to lattice sites in the neighborhood 
of each particle as a function of the particle position. 
Suppose that each charge q a at continuous position r a is 
spread over lattice sites {R^}. The charge on site Rfe is 
given by g a s(r a ,Rfc). We call s the charge assignment 
function which has the property 

^s(r ,Rj-) = 1. 

3 

Many choices of assignment function are possible and the 
choice may greatly affect the performance of the algo- 
rithm. Clearly, as the number of lattice sites over which 
each charge is distributed increases, the computational 



cost increases also. As for the accuracy, the overlap of 
lattice based charge clouds results in a deviation from the 
continuum potential for charge separations on the scale of 
the charge cloud. Because, for high accuracy, one should 
correct the energy for charges with overlapping charge 
clouds, there is a greater computational cost for higher 
order assignment functions. For Fourier transform meth- 
ods, including most molecular dynamics methods Q, ac- 
curacy can increase as the radius of the charge cloud 
increases. This is because narrow charge clouds have a 
large spread in Fourier space and higher frequency modes 
can act as aliases for lower frequency ones in the sense 
that modes of different frequency can appear the same on 
the scale of the lattice if their values at the lattice sites 
coincide. To avoid this aliasing error, it is better, for 
these simulation methods, to have larger charge clouds 
so that frequency aliases have a smaller amplitude. We 
use a third order scheme due to Hockney and Eastwood 
[To| that spreads the charges over nine vertices in two 
dimensions. If Dx and Dy are the distances in units of 
the lattice spacing a of the charge q to the vertex wo,o 
closest to it, the weights for charge assignment in the x 
direction are 

Wl x = {Dx - 0.5) 2 /2 

Wfi = 0.75 - Dx 2 (6) 
VF+i = (Dx + 0.5) 2 /2 

with similar formulas for the y direction assignment. The 
charge assignment function s is a product of pairs of 
weights as shown in Figure [TJ 

The total energy of the system in terms of the Green's 
function is 

U = \ E E E 9 2G ( R < - Ki)a(r a ,R i )a(T b ,B. j ). (7) 

r a r b i.j 

Alternatively, it is the sum of all the squared link electric 
fields 

2 N x -1 N v -1 

U =Y E E (E 2 ((n x ,n y ),x)+E 2 ((n x ,n y ),y)). 
n x =0 n y =0 

(8) 

This interaction is the one generated by the local sim- 
ulation algorithm. At charge separations greater than 
about one lattice spacing (depending on the charge as- 
signment scheme used) this interaction energy is an ex- 
cellent approximation to the Ewald potential. At short 
range the interaction is weaker than in the continuum 
problem. In fact the limit of the interaction energy as the 
separation of two charges tends to zero is finite. For accu- 
rate simulation this must be corrected. For each charge 
move, charges in the neighborhood of the charge that is 
being moved are identified using the linked list method 
of Hockney and Eastwood [ioj . The discretized inter- 
action energy from equations {7} (with the sum suitably 
restricted) and (O is subtracted off for each pair of neigh- 
boring charges and the Lekner potential V from equation 
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([3]) or a good short range approximation to the Ewald po- 
tential is added on in its place. A suitable approximation 
for small charge separations r is 



Usrc = ~7T lnr 

Z7T 



2 2 

~4A 



(9) 



where A = L x L y , which is derived in an appendix and 
presented with a calculation of the deviation from the 
Lekner potential for different charge separations. For 
very high accuracy or for low temperatures, the Lekner 
interaction is more appropriate. The linked list method 
is a local method so that the scaling of the local algo- 
rithm is preserved. The periodic cell is divided up into 
chaining cells - regions of 2 x 2 plaquettes in our simu- 
lations. At the beginning of the simulation, the charges 
are given chaining cell labels and this information is put 
into a linked list. The linked list is composed of two ar- 
rays; the first gives a charge label given a chaining cell 
label - in other words it gives the label of one charge in 
a given chaining cell; the second array has particle labels 
as its contents and also as its indices in such a way that 
given a particle label, another particle in the same chain- 
ing cell is found using the previous element as an index, 
the new element being zero if there are no more charges 
in the chaining cell. For a more detailed explanation of 
this computer algorithm we refer to the book by Hockncy 
and Eastwood [ldj . 

There is another problem associated with working on a 
lattice: the interaction includes finite self-energies for the 
charges: because each point charge is spread over several 
lattice sites, it interacts with itself with finite energy. 
The self-energy in this context is a one particle potential 
energy that has regular wells at positions relative to the 
lattice sites. This energy must be subtracted off, once 
again using (JT]) (with the sum taken over the sites over 
which each charge is spread) and ((5]). This subtraction is 
essential; without it, one finds that the charges become 
trapped within the lattice-induced potential wells at low 
temperatures, while higher temperature results are also 
strongly affected by this potential. 



III. LOCAL METHOD 
A. Continuum update 

The novelty of the local algorithm is in the way the 
field configuration changes. When a charge moves, the 
field changes in the vicinity of that charge so that Gauss' 
law is preserved. We first consider charge moves in the 
continuum limit to illustrate some features of the dy- 
namics. Because, in this subsection, we are not directly 
concerned with computations, we work with SI units and 
in three dimensions. 

The electric field at time t due to particle motion can 
be written as 

E(r,t) = E(r,to)- V- [ dnS{r - n(t)) 

A C J 



where E(r, to) is the purely longitudinal electric field at 
time to when the charges are in their starting positions. 
The electric field, according to this equation, evidently 
varies locally as the charges move. Each charge leaves a 
trail along its path as pointed out by Levrel and Maggs 
(TH . If we take any closed surface around some charges 
at to, Gauss' law is satisfied by the definition of E(r, to). 
As time increases, charges may enter or leave this closed 
surface, whereupon the delta function ensures that flux is 
chopped out or added in at the cutting point to preserve 
Gauss' law. If F is the closed surface 



■1 <f dS- [ dr t 6(r-Ti(t)) 
£o Jf J 



iq/eo if the charge cuts F an odd number of times 
otherwise 



where the sign on the right-hand-side depends on whether 
the charge is originally inside (minus sign) or outside F. 
The electric field trail term can be written as 



dr,V 2 



f 



4tt€ J ~" Vl r - r i(*)l 

in three dimensions. Using the identity 

V 2 V = V(V • V) - V x V x V 

for vector field V we notice that, in accordance with the 
Helmholtz theorem, the electric field can be split into 
longitudinal and transverse parts 



E = V0 + V x Q 



where 



AlTEo 



-V • / dr 



1 



#{*(*<>)}) 



and 



47re 



-V x dr 



|r-r*(i)| 



We can go further and write the scalar and vector poten- 
tials as 

, M ST q f dri " ( r ~ r *( f )) 

^)=Li^i |r-r,(t)|» 



£ 



47T60 



dri • Vi 



r - n 



mn(t )}) 



which is recognizable is the ordinary electric scalar po- 
tential while 



Q = £ 



q f d^ x (r - ri(t)) 



Atteo 



|r-ri(t)|3 



which is the Biot-Savart law. So, the electric field at each 
time t is the Coulomb field due to the charges at time t 
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positions, and a circulation V x Q where Q is equivalent 
to the magnetic field produced by current-carrying wires 
along the trails left by the charges as they move - the 
current being qj 

We notice that if a charge makes a closed loop, there 
is no change in the longitudinal field as we should ex- 
pect. Secondly, the transverse field is always zero for the 
closed loop except on the loop itself, in accordance with 
V x Q being proportional to the current density in or- 
dinary electromagnetism. When the path is not closed, 
this relationship breaks down; in terms of the current- 
carrying wires analogy, the cndpoints imply that charge is 
not conserved so there is an additional magnetic field con- 
tribution, the curl of which is identical to the Coulomb 
field of the charge at the endpoints. We have demon- 
strated, that pure transverse field updates occur entirely 
by the motion of charges in closed loops. 



Dx 



Dy 



B. Off-lattice charge move update 



B 



The O(N) Monte Carlo charge update requires us to 
maintain the discretized version of Gauss' law, equation 
((4]) , at each charge move by changing a few links as pos- 
sible for a given charge assignment function. The energy 
difference computation at each step from equation ([8]) in- 
volves only the link fields that have changed. There are 
many possible ways of updating the electric field locally. 
When charges are completely confined to the lattice sites, 
then hopping of charge q across a link requires us to alter 
the electric field E on that link by E — » E — q/a to pre- 
serve the flux constraint. The charge leaves a trail of elec- 
tric field, as we have already seen in the continuum case. 
I have used the following field update procedure for off- 
lattice simulations ij. Given the restrictions on charge 
motion and the assignment function, the number of ver- 
tices affected by the charge move must be either nine or 
twelve. Suppose that the nearest vertex to the charge Ro 
is the same before and after the move so that only nine 
vertices are affected; the other case is identical in princi- 
ple. For each affected vertex R(p, q) = Ro +px + gy with 
p,q = —1,0, 1, the change in the charge is Aq(R(p, q)). 
Beginning at the vertex R(— 1,— 1), we trace out an 5* 
shape, updating the links along the path so that, at each 
site along the path, Gauss' law is satisfied. Figure Q] 
shows the general idea. 



C. Other features 

In continuum terms, the local update is responsible for 
generating an electric field of the form 

E = -V(f> + V x Q + E fc (10) 

with Q and uniform field non- vanishing at finite tem- 
perature. By definition, the potential produces an elec- 
tric field with no net background field E&, but the local 



FIG. 1: Charge update for charge move in neighborhood of 
a single plaquette. Completed link updates are shown in red. 
Those to be updated in green. The original charge posi- 
tion from the nearest vertex is given by Dx and Dy and the 
charge assignments on vertices 1,2 and 3 are initially given 
by Wl^W-x, W\W§ and W^Wf respectively where the 
charge assignment weights W are given in equation ©. If, 
when the charge moves, the change on vertex i is Aq(i) then 
the link field denoted A is updated by Ea — * Ea + (Aq(l)/a) 
and for link B E B -» E B + ((Ag(l) + Ag(2))/a) and so on. 



field update gives rise to fluctuations in this field. To 
see this, consider a charge that hops from one lattice site 
to a neighboring site. The simplest local field update 
produces a field trail that would restore the charge to 
its original position at zero temperature. If the uniform 
background were zero before the charge hopped between 
sites, then because the background field is the sum over 
the vectors at each site divided by the total number of 
links there must such a field in the final state. This field 
can be spread over the lattice by further field updates. 
From this kind of argument, we see that the background 
field that is generated depends on the distance of the 
charges from their original positions. More precisely, at 
time t, a single charge will be responsible for a uniform 
electric field that is proportional to q(r(Q) — r(t)) where 
r(0) is the original position of the charge q and r(t) is its 
position at time t. It should be emphasized that r is mea- 
sured without any regard for the periodicity of the cell. In 
other words, the electric field increases as a charge winds 
around the cell so the charge experiences a restoring force 
towards its original position. With many charges, the 
uniform field is proportional to X)ift( r i(0) — r iW)- So 
we should expect that fluctuations in the uniform field 
will cause charges to lose sense of their original positions. 

It has been argued that the local update can produce, 
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even after averaging, an interaction of the more general 
type considered by Leeuw and Perram - a periodic po- 
tential and a term that depends on the dipole moment of 
the cell Q and which is equivalent to the presence of a 
uniform electric field within the cell. We recall from Sec- 
tion [II] that the periodic interaction is recovered from a 
spherical cluster of identical cells by cladding the cluster 
with a perfect conductor. If the cladding has finite dielec- 
tric constant e, a dipole term is introduced as in equation 
fl2J). The dipole term in the cell is bounded because the 
cell is finite. In contrast, the uniform field produced by 
the local is unbounded and does not depend on the abso- 
lute dipole moment unless the dipole moment of the cell 
is zero at the beginning of the simulation. The argument 
is that at high temperatures, charge winding causes 
the uniform electric field to average out during the simu- 
lation. At low temperatures though, the field does make 
a difference to the effective interaction after thermally 
averaging. 

One can see that the uniform field makes no difference 
to thermal averages provided it is sampled independently 
of charge positions. So one way of ensuring infinite e is to 
introduce an independent background field update. This 
can be implemented locally. It has the disadvantage of 
slowing the simulation down. A more efficient means is to 
keep track of the uniform field locally and to subtract off 
its contribution to the energy. This can be done because 
the uniform field energy decouples from other contribu- 
tions even before thermal averaging. We shall investigate 
the influence of the uniform field further in the next sec- 
tion. 

Although charge updates do sample transverse electric 
field configurations, we expect to have to sample these 
circulations independently of the charge moves. For the 
circulation sampling, a plaquette is chosen at random 
from a uniform distribution over the whole lattice. A 
number A is then selected from a uniform distribution 
over the range [— v, v\ where v is chosen so that the move 
has a 50% acceptance rate. The circulation of the electric 
field around the chosen plaquette is changed by altering 
each of the links of the plaquette by A as shown in figure 
[H The energy difference calculation therefore involves 
only four link variables. 



D. Literature 

Monte Carlo with O(N) scaling by the introduction 
of transverse field degrees of freedom may have charges 
confined to the lattice as introduced in [l|,[ll(, or allowed 
to move off-lattice with some smooth charge assignment 
onto the lattice vertices Q which we have followed in our 
simulations. 

At first sight, we would expect that the transverse field 
degrees of freedom would have to be allowed to relax 
across the whole system for each particle move. Although 
the scaling of the speed varies as order N when only parti- 
cle moves are considered, one would naively expect order 
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FIG. 2: Plaquette circulation update at lattice site (ij). The 
electric field divergence at each site is preserved. 

N 2 scaling when transverse field sampling is included. It 
turns out, however, [12] that one can sample the trans- 
verse fields much less frequently than one would expect 
without spoiling the results. This is a surprising result. 
In connection with this issue, Maggs and coworkers (l3j 
have drawn attention to the fact that the transverse elec- 
tric field is sampled by charge moves without requiring 
independent plaquette updates. In this paper we expand 
on this observation. 

It has been found that for charges confined to a lat- 
tice, charges tend to freeze in their original positions. 
Two papers are devoted to various ways of removing this 
obstacle in on-lattice simulations [ll|, [l4| by improving 
the efficiency of circulation updates or by spreading the 
charges over several sites. For off-lattice simulations, this 
problem can be present but one can avoid it by taking 
the step size to be small or by spreading charges over 
more sites. 

The local update method has been developed by Maggs 
and collaborators for molecular dynamics [15l |l6[ and in 
work by Pasichnyk and Diinweg [16j |. the competitive- 
ness of local molecular dynamics was remarked on in a 
comparison with a global O(N) method. 

The method has been applied in various test cases in 
the original papers. More recently however, it has been 
used by Maggs to carry out Monte Carlo simulations of 
a model of water [l7l |. 

IV. RESULTS OF COMPUTER EXPERIMENTS 

A. The one-component plasma in two-dimensions. 

Our simulations focus on the one-component plasma 
in two dimensions. This is a system of particles, all with 
charge q, on a fixed neutralizing background. In two 
dimensions the interaction between the charges is loga- 
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rithmic. The phase diagram of this plasma depends only 
on the coupling T where 2irT = q 2 fi. For T = 2 the model 
can be solved exactly for the free energy and the n-point 
correlation functions. At around T = 140, Monte Carlo 
simulations using periodic boundary conditions have re- 
vealed the existence of a first-order transition from a liq- 
uid to a crystalline phase. Most of the numerical work 
in the following focuses on this model. Firstly, the ex- 
act solution allows us to perform an objective test of 
the accuracy of the lattice based methods. Secondly, the 
plasma has long-ranged correlations for T > 2 which form 
the basis for a stronger test of the achievable accuracies. 
Thirdly, the model has a simple formulation and is suit- 
able for making clean tests of other properties of the new 
method including the effect of varying the rate of inde- 
pendent circulation updates. 



B. Simulations of the 2D OCP 

I have carried out Monte Carlo simulations of the one- 
component plasma in two dimensions using the local elec- 
trostatic algorithm discussed at length in the previous 
section and also a continuum method using the Lekner 
form for the pairwise interaction. The simulation cell for 
the discrete methods is divided into 32 x 32 plaquettes 
and each chaining cell is a 2 x 2 plaquette region. The 
short range correction algorithm extends over a 6 x 6 
plaquette square with the charge that is moved in the 
central chaining cell. I used the approximation to the 
Lekner potential ^ as the correcting potential. All the 
results in this subsection were taken with 120 charges in 
the simulation cell. The lattice spacing a was chosen so 
that the number density of the particles was unity. The 
rate of plaquette updates was usually taken to be 40 for 
each charge move. Differences between the results ob- 
tained here and the correct thermal averages arise from 
small differences in the potential energy for small charge 
separations because an approximation to the correct po- 
tential is used (see appendix [A| and also from imper- 
fect sampling when the relaxation times are large which 
happens towards the lower temperature range we have 
investigated. If the approximation to the exact interac- 
tion were insufficiently accurate, one could use the exact 
potential instead. 

For various values of the coupling T, I have computed 
the pair correlation function g(ji, r^) which is defined as 

ff(ri,r 2 ) 

where Pjv(ri, •••> r Jv) is the thermal probability distri- 
bution function for the 2D OCP with N harmonically 
confined charges with the number density set to one. 
This quantity is proportional to the probability of find- 
ing charges in the neighborhood of positions ri and r 2 . 



Because the liquid, in the thermodynamic limit, is ho- 
mogeneous and isotropic, the pair correlation function 
depends only on r = |ri — ra|. This definition for the 
pair correlation function has the property that it tends 
to unity for large values of its argument. For T = 2, 
Ginibre [HI has found g{r) exactly to be 

g(r) — 1 - exp(-7rr 2 ). (11) 

This exact result exists because the partition function 
for N charges is identical to the normalization of a wave- 
function of free fermions with Gaussian wavepackets. In 
Figure [31 the exact Gaussian curve is plotted with results 
from the Monte Carlo simulations. There is only weak cor- 
relation for small charge separations; otherwise g(r) is 
featureless. Both methods reproduce the correct behav- 
ior extremely well; for the case of the lattice method, 
the correlation function is correct at small separations 
between the charges provided the short range correc- 
tion method is used. Without correcting the potential 
for small charge separations, the radial density does not 
tend to zero as r tends to zero, but is close to the exact 
solution for charge separations greater than about one 
particle spacing. 

At higher values of T, we make a direct comparison of 
pair correlation functions from Monte Carlo data. There 
is excellent agreement between the Monte Carlo meth- 
ods for r < 100 as shown in Figured) As we should ex- 
pect, correlations become more pronounced and longer- 
ranged as r increases. For lower temperatures, discrepan- 
cies between the radial densities for simulations using the 
Lekner potential and the local simulation appear. They 
are small at T = 100 but become very pronounced for 
T ~ 140 which roughly corresponds to the reported on- 
set of the crystalline phase. Correlations are consistently 
weaker for the local method. 

By increasing the range of the short-range correction, 
the results for the local method can be brought arbitrarily 
close to those for the Lekner summation at the expense 
of slowing down the energy computations. For example, 
if we extend the range of the short-range correction from 
6 x 6 to 12 x 12 plaquettes, we can obtain a perfect match 
between the Lekner correlation function and the local 
method result as shown in Figure [5] 

Because the agreement between the Lekner potential 
and local method results is excellent even down to very 
low temperatures, we rule out the possibility that the 
background field can make a significant contribution to 
the potential. As a further test, we have included an 
explicit uniform field sampling term in our simulations. 
The results are shown in Figure [51 The pair correlation 
functions match the Lekner result with or without the 
extra update when the 12x12 correction region is used. 
A difference does show up within the smaller correction 
region but the result with the extra update is less corre- 
lated than the result without the field update. I believe 
that one cannot rule out the winding of the electric field 
at any temperature provided the charges are mobile and 
that therefore, the uniform field will average out in most 
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FIG. 3: Pair correlation functions for the 2D one-component 
plasma at F — 2 with (red) and without (green) short-range 
potential correction. When the short-range potential is cor- 
rected the exact solution (also plotted) and simulation result 
overlap. 




FIG. 4: Pair correlation function g(r) for the 2D one- 
component plasma for T equal to 20,40,60,80,100. Results 
plotted for simulations with local update, and with the po- 
tential energy computed as a Lekner summation. 



cases. 

I have also computed the average energies of the 
charged particle system using the Lekner summation as 
a further test of the local method. Of course, global en- 
ergy calculations at each step defeat the object of the 
local algorithm, which is most efficient if we are only 
interested in charge correlations. Anyway, the energies 
from this work agree temperature for temperature with 
results obtained by Leeuw and Perram Q to their degree 
of accuracy. 



C. Autocorrelation times of charges 

In this section and the following, we present various 
autocorrelation times denoted r which are characteristic 



FIG. 5: Pair correlation function g(r) for the 2D one- 
component plasma for F = 140. Five curves are shown. The 
gray curve is for a global continuum method. The green and 
blue curves that match it were found using the local method 
with a short-range correction extending about 1/3 of the to- 
tal distance across the system from each charge. They are 
distinguished by either including independent uniform field 
sampling (green) or not (blue). The two less correlated curves 
(red and black) are respectively with and without uniform 
field sampling. They are less correlated than the other three 
curves because the corrections to the potential extend only 
half of the distance from each charge compared to the more 
correlated case. 



times of the autocorrelation function 



A(k) = 



(O t O l+k ) - {Oif 

{oj) - <o>) 2 



where 0, for i = 1, . . . , N is a sequence of data obtained 
from the Monte Carlo simulation and the angled brackets 
denote a thermal average. These times provide a mea- 
sure for the rate of relaxation of different observables 
over the course of the simulation. They depend on the 
update method, on the nature of the particles and their 
interaction and on the choice of O. We measure these us- 
ing both a binning procedure p^ . [20| and by integrating 
the autocorrelation function. We find the numbers are 
consistent with one another, though the former method 
is preferable because it signals convergence to r. We 
present results in this section for two typical observables 
at various temperatures. After an equilibration process 
from a hot start, the algorithm is allowed to adjust the 
maximum step size Ar that each charge can make from 
its original position so that the acceptance rate for charge 
moves is 50% [251 ] . The sampling is performed with this 
fixed maximum step size. A similar and simultaneous 
procedure is adopted for the plaquette update step in the 
local simulation. Because we are concerned with charge 
observables in this section, we measure autocorrelation 
times in units of charge move trials; in other words inde- 
pendent transverse field updates are not included in the 
estimates. 

Figure [5] shows a rise in the energy autocorrelation 
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FIG. 6: Autocorrelation times for the energy on a logarith- 
mic scale varying with F. Fhe circles are data for the local 
method, squares for the Fekner simulation. 



time as T increases; as the temperature is lowered, en- 
ergy fluctuations become more improbable so the step 
size Ar falls to retain the fixed acceptance rate. The 
autocorrelation time falls coincidently. Over the range 
r = 2 to r = 80 there is roughly a tenfold increase in 
r for both the Lekner and local methods. The autocor- 
relation times are consistently much larger for the local 
method than for the Lekner method. That the local sim- 
ulation Ar should be smaller than that for the Lekner 
simulation is to be expected on the following grounds. In 
the absence of charges, a single charge has to overcome 
an energy barrier, on average, to make a move using the 
local method. This is because it has to move in its own 
field which relaxes locally to maintain Gauss' law. This 
can be seen most easily in the case of a single charge 
confined to a lattice vertex. If we ignore field circula- 
tion changes which only produce fluctuations about the 
average that we consider here, a step in the position of 
the charge takes a single link field from E to E — (q/a) 
where E = q/4a (in 2D) on average; the energy differ- 
ence is then q 2 /A > 0. In contrast, a single charge in a 
periodic cell in the Lekner simulation has no barrier to 
overcome to make a move on average. The qualitative 
observation that the energy barrier for a given step size 
is higher in the local method than in the global method 
is true regardless of the number of charges present in the 
cell. 

I have also computed the structure factor for k vec- 
tor close to the first reciprocal lattice vector for our two 
different simulation methods. As correlation develops in 
the sample, a peak forms in the structure factor at recip- 
rocal lattice vectors. The structure factor at these values 
of k are the slowest to relax so we have concentrated on 
high temperatures compared to the reported phase tran- 
sition temperature to ensure good results. Once again, 
we find that r is much larger for the local simulation 
method than for the global method as can be seen from 
Figure [71 which uses a logarithmic scale to represent the 
autocorrelation times. 



FIG. 7: Variation of log 10 (r) of the structure factor for k vec- 
tor (10, 6) (close to the first reciprocal lattice vector) with in- 
verse temperature V in the liquid phase which exhibits short- 
range correlations. The local simulation (circles) relaxes more 
slowly than the global simulation method with a Lekner po- 
tential. 

D. The relative speed of the local algorithms 

In previous works, the O(N) scaling of the local 
method has been pointed out. But, such scaling can- 
not be of much use if the scaling prefactor or, in other 
words, the absolute times are prohibitively large even for 
small systems. In this section, we shall see that in com- 
parison with the Lekner summation simulations, the local 
method fares very well even when the differing autocorre- 
lation times are included. The autocorrelation time pro- 
vides a measure for the efficiency with which simulation 
methods can sample configuration space. Specifically, if 
a single run gives a sequence of n values of some observ- 
able with variance Oq. , the variance of the data about the 
ensemble average for that observable is given by [lj| [2(| 



We use the autocorrelation time of some observable to 
find the number of sweeps required to bring the variance 
about the ensemble mean down to some value and hence 
obtain the absolute simulation times taking into account 
the differing charge relaxation times. Different observ- 
ables have different autocorrelation times in general but 
to make a comparison of absolute speeds in this way we 
have to choose one. We have chosen 5(k = G), because, 
as a charge correlation function, it is a quantity that is 
both interesting and which the local algorithm is partic- 
ularly suited to computing. Also, the reciprocal lattice 
vector is likely to give the slowest relaxation of charge 
observables, so the test will provide something of an up- 
per bound on times. It should be emphasized that the 
autocorrelation time for sufficiently high temperatures or 
sufficiently large system sizes, is independent of system 
size provided it is measured in Monte Carlo sweeps. So, 
the comparison made in this section does not affect the 
scaling of the algorithms. This is not the case close to 
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phase transitions, which would require a separate analy- 
sis, or for very small system sizes. 

I have measured the time taken for both the local and 
Lekner algorithms to complete a given number of particle 
moves for two system sizes N — 120 and N = 1000. This 
time depends, of course, on the code, the compiler and 
the machine that is used so we don't expect the numbers 
presented here to demonstrate any more than trends in 
the absolute and relative speeds and rough estimates for 
the absolute speed. The Lekner method scales as 0(N 2 ) 
but it has the advantage that it is easy to implement and 
the time taken for a charge move should be smaller than 
in O(N) molecular dynamics methods for typical system 
sizes. 

The code with local updates has a short-range poten- 
tial correction routine with a 6 x 6 plaquette region within 
which the correction is performed. The simulation cell is 
a 32 x 32 mesh. The number of plaquette updates per 
particle move is maintained at 40 so that the mobility is 
close to maximum as we shall see in Section IIVEI The 
cell size is chosen so that the charge density is one. 

The raw speed of the local code before including the 
correlation in the data is obviously much greater than 
that of the global code because the overheads are so much 
smaller when just a few links need to be changed. We 
also find that the times scale as N for the Lekner code as 
expected. The time interval for a fixed number of steps in 
the local code does vary with N only because the number 
of particles captured by the short-range correction rou- 
tine increases with N. In fact, the short-range correction 
makes a very significant contribution to the speed: for a 
cell with 200 charges, switching off the correction routine 
halves the time interval. When the number of plaquettes 
increases holding the number of particles per plaquette 
constant, the speed of the local update code does not 
change significantly. 

In the previous section, the variation of r for S(k = 
G) with temperature was given. We have also measured 
the mean and variance of the structure factor for each 
temperature from the data. Suppose we want to estimate 
S(k = G) to within 1%. We obtain the required number 
of configurations n from equation (|12[) . 

Table U gives, for seven temperatures, the autocorrela- 
tion time of the first reciprocal lattice peak of the struc- 
ture factor (which has already appeared in Figure [7]), the 
variance in the configuration values about the mean and 
hence the number of sweeps required to get the error of 
S(k = G) to 1% of the mean value. From our measure- 
ments of the speed of each piece of code, we estimate 
the minimum amount of time the global and local sim- 
ulations should be allowed to run to obtain data of this 
quality for N = 120 and 1000. The effect of including 
the autocorrelation times is contained in the two lowest 
rows. The local code remains the faster of the two, for 
both system sizes and for all temperatures but inclusion 
of the relaxation times has narrowed the gap consider- 
ably. For r = 20 and N = 120 the necessary times are 
brief (as we would expect at these fairly high tempera- 
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7.7 
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11.0 
80.7 


17.8 
162.6 


23.7 
199.4 


27.8 
325.6 


42.8 
446.0 


53.9 
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Simulation time 
120 charges/mm 


10.7 

4.9 


15.4 
8.6 


24.9 
17.3 


33.2 
21.2 


38.9 
34.7 


60.0 
47.5 


75.0 
68.3 


Simulation time 
1000 charges/hr 


12.6 
3.2 


18.1 
5.7 


29.2 
11.5 


39.0 
14.0 


45.7 
23.0 


70.5 
31.5 


88.7 
45.3 



TABLE I: Table of speeds of the Lekner and local simulation 
methods taking into account the fact that their autocorrela- 
tion times differ. We used a 1004MHz, AMD Athlon x86 64 
bit processor with Linux operating system and suitably opti- 
mized gec compilation. The upper figure in each box is for 
the Lekner summation Monte Carlo, the lower number for the 
local method. 

tures compared to the phase transition) and similar. The 
scaling of the speeds with the number of charges implies 
a larger difference for N = 1000 but the absolute times 
are becoming inconveniently large. 

E. Transverse field relaxation 

The unique feature of the local simulation method is 
the introduction of transverse field degrees of freedom. 
These degrees of freedom should scale roughly with N, 
the number of charges in the system to preserve the ac- 
curacy of results. One might expect that the transverse 
electric field should be allowed to relax all the way across 
the system at each charge move, with the result that the 
algorithm should scale as 0(N 2 ). In practice, we find 
that a very low rate of independent plaquette updates 
O(l) is all that is required per charge move to obtain 
good sampling at intermediate temperatures 2 < T < 140 
and, at higher temperatures, plaquette updates can be 
switched off altogether. At lower temperatures, we find 
that one must sample the transverse field at a higher rate 
per charge move. 

In this section, we investigate these properties. To 
do this, we have measured the autocorrelation time r of 
the background uniform electric field as the rate of pla- 
quette updates is changed. The uniform field does not 
depend directly on the plaquette circulation; it varies as 
charges move by the local update method as discussed 
in Section [TlJ so this quantity provides a measure for 
the charge relaxation times that is particularly straight- 
forward to investigate compared to relaxation times of 
charge correlation functions. 

Figure[H]shows how r varies when the plaquette update 
rate is varied. The horizontal axis in the figure is 

Rate of plaquette updates 

= a. 

Rate of particle updates 

We find that, as a increases from zero for intermediate 
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temperatures, r decreases sharply and reaches a vale. 
Over the temperature range 5 < T < 140 that we have 
investigated, the behavior of r is so weakly dependent 
on temperature that the data almost collapse. A similar 
behavior is seen in the mobility of charges, (defined here 
as the gradient taken from the plot of RMS distances of 
charges from their original positions against simulation 
time measured in Monte Carlo sweeps). The mobility is 
close to maximum for a very low rate of plaquette up- 
dates: just one update per particle move on average al- 
though there are 32 x 32 plaquettes and roughly one par- 
ticle for nine plaquettes. We have also looked directly at 
the autocorrelation times of different modes of the trans- 
verse field with the same temperature independence even 
close to the transition. At higher temperatures T < 20, 
we see a significant change in the behavior of r for the 
uniform field. The same figure [8] shows data for T = 2 
for which r falls below the lower temperature data at 
r = 5 at small a. We have also been able to obtain the 
a = point for T = 2 which is too large at lower tem- 
peratures to obtain accurately. This high temperature 
behavior matches the qualitative result that good aver- 
ages are obtained at high temperatures T < 2 without 
plaquette updates. 

We can understand this behavior qualitatively as fol- 
lows. We have already pointed out that in the absence 
of plaquette updates, the electric field is altered by the 
motion of charges, which leave trails of altered field along 
their paths. The effect of the trail on each charge is to 
restore it to its original position. So, at least, at low 
temperatures, we expect charges to be trapped in their 
original positions with a correspondingly long autocorre- 
lation time of the uniform field which oscillates with the 
charges. When the electric field is allowed to relax by 
changing the circulation around plaquettes, the string is 
smeared across the system and the restoring force on each 
charge falls on the average and hence the charges become 
mobile. There must be some limiting mobility at a fixed 
temperature. We expect that smearing of the trails can 
also take place purely by the motion of charges - an effect 
which is presumably enhanced by spreading charges over 
more lattice sites. But this smearing of trails by charges 
is contingent on at least some independent plaquette up- 
dates. At high temperatures, charges can overcome the 
potential barriers around their initial positions so r for 
the uniform field is small. 

We would like to understand why the maximum mobil- 
ity is attained with such a low update rate at intermedi- 
ate temperatures. We shall see that a simple probability 
model for the update process exhibits a similar feature. 
Suppose there are N plaquettes, each one inhabited, on 
average, by a charge. At each step, the probability of 
choosing a plaquette is P p and the probability of choos- 
ing a charge P c . The ratio P p /P c = ct. So 



P, 



1 



Pc 



1 



can move freely (it does not leave an electric field trail). 
Otherwise, it is trapped. The number of free plaquettes 
is N F . The following processes are possible 

1. Plaquette move frees plaquette: N F (ti_i) + 1 = 
N F (ti) 

2. Plaquette move traps plaquette: Nf(U-i) — 1 = 
N F (ti) 

3. Particle move on trapped plaquette: N F (ti_i) = 
N F {t t ) 

4. Particle move on free plaquette: N F (ti-i) — 1 = 
N F (U). 

The fourth rule that a charge moving on a free plaquette 
causes the plaquette to become trapped has its basis in 
the fact that a moving charge leaves a line of electric field 
when it moves. The probability distribution function of 
the number of free plaquettes satisfies a master equation 



dP(N F ,t) 
dt 



J {N F ) — J (N F + 1, i) 



where 



1 



A plaquette is said to be free if a charge on that plaquette 



q, jv - Np + 1 N F 
J(N F ) = — jf^P(N F - 1, t) - -^P(N F , t). 

We solve for the stationary (time-independent) probabil- 
ity distribution in the usual way 21]. Form the summa- 
tion 

M-l 

Y, [J(i + l)-J(i)} = J(M)-J(0). 



The boundary conditions on the rates give J(0) = and 
hence 

a N - Np + 1 

P(N F ) = — - ^^—PiNF 1). 

a + 1 Np 

We compute the mean number of free plaquettes (N F ) as 
the ratio a is changed. The mean number of free plaque- 
ttes is essentially equivalent to the number of free charges 
and is therefore indicative of the mobility. We find that 
(N F ) tends rapidly to a maximum of half of TV as a is 
increased (see Figure [3]). Around 10 plaquette updates 
per particle move are sufficient to bring (Np) within 5% 
of maximum just as is the case for mobilities in the simu- 
lations presented above. Of course, this model is hugely 
oversimplified. There are many ways in which the sim- 
ulation updates differ from the processes in this model. 
For example, the simulation uses a Monte Carlo test to 
carry out updates, charges are spread over four plaque- 
ttes and we have simulated at a particle density of about 
1 particle to 9 plaquettes (neglecting spreading). Also, a 
charge can erase a trail by moving along it: a trapped pla- 
quette can be freed up by charge motion. Nevertheless, 
the model captures the important feature of the changing 
mobility: if more than one plaquette is freed per particle 
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FIG. 8: Variation of the integrated autocorrelation time of 
the uniform background field with plaquette update rate mea- 
sured in average number of plaquette moves per particle move. 
Results for four different temperatures are shown: V = 2,5, 
20,40, 60 and 140 (circles,stars,crosses,squares, diamonds and 
triangles respectively). The definition of a is given in the 
main text. The inset shows the gradient of RMS distance of 
particles from their original positions varying with V. 
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FIG. 9: (Nf) as a function of a showing the rapid approach 
to 512 where there are 1024 plaquettes in total. 



move, Nf increases until limited by the unfreeing oper- 
ation. We therefore expect the limiting number of free 
plaquettes to be half of the total. 

We turn now to the behavior at low temperatures 
around the transition temperature. We find that a high 
rate of plaquette updates is necessary to achieve the same 
degree of correlation one finds with the Lekner summa- 
tion method. For T — 120 and 140, I have obtained pah- 
correlation functions at different values of a in the range 
a = 2 to a = 100. All the T — 120 results are identical 
and match the pair correlation function using the Lekner 
Monte Carlo method. The T = 140 results differ, with 
convergence onto the Lekner result for around a = 40. I 
have not investigated this further but it is worth point- 
ing out that the slower convergence to the correct result 
as a is varied happens close to the transition tempera- 
ture where one would expect charge relaxation times to 



be longer than in the fluid phase. There is no remark- 
able behavior close to the transition temperature in the 
uniform field relaxation time nor in the transverse field 
relaxation time. For practical purposes at such low tem- 
peratures, the unbiased plaquettes updates arc inefficient 
and it may be sensible to use a worm update to sample 
transverse degrees of freedom [13j which involves the cre- 
ation of a pair of equal and opposite charges, at intervals 
during the simulation, which wander around the simula- 
tion cell and eventually annihilate. 

We conclude this section with three observations. The 
first concerns the relationship between the mobility (de- 
fined above) and the autocorrelation time of the uniform 
field t. We have found that as t falls, the mobility rises. 
For low a, this is consistent with the picture of charges 
trapped by their electric field trails. However, the mobil- 
ity is also temperature dependent as can be seen in the 
inset plot of Figure [5J At low temperatures, the mobil- 
ity is small although the uniform field may relax quickly. 
So, r appears to be independent of the mobility. Notice 
however, that the greatest temperature change in the mo- 
bility and in r occurs for T < 20 whereas for T > 20 both 
are only weakly temperature dependent, so they are more 
(anti) correlated for low a. 

The second observation is that the curve of uniform 
field r(a) depends on the number of charges N and 
the number of plaquettes P only in the ratio P/N. If 
the number of plaquettes increases, holding the number 
of charges constant, the autocorrelation time increases 
because the fraction of plaquettes updated per particle 
move drops. When N increases in proportion to P, the 
data collapses onto the curve in Figure [H 

Finally, the results described in this section lead to the 
conclusion that if charges can overcome the impetus to 
retrace the electric field trails they produce (which will 
typically be at high temperatures), then thermal aver- 
ages will be electrostatic averages even in the absence of 
independent plaquette updates. After all, the Metropolis 
algorithm guarantees that averages from computer gen- 
erated configurations will converge to thermal averages 
(which in this case are electrostatic averages) provided 
that all of configuration space can be sampled in princi- 
ple. In practice, reasonable mobilities can be achieved at 
high temperature with no plaquette updates. At lower 
temperatures, a few plaquette updates are required to 
overcome bottlenecks while electric field trails are altered 
by charge motion and charges can cooperate in freeing 
other charges provided they are dense enough that charge 
clouds overlap to some extent. 

Although charge dynamics depends on electric field cir- 
culation and vice versa, the factorization of transverse 
and longitudinal energies ensures that these contribu- 
tions to the fields are importance sampled independently 
and simultaneously. Relaxation of transverse fields across 
the periodic sample at each charge step is unnecessary 
because we are interested in averages over many config- 
urations. It is sufficient that relaxation can occur over a 
number of charge steps. At the end of the simulation with 
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no information discarded, one can 'project' onto each link 
transverse field to find it has been sampled adequately. 
Each transverse field is a harmonic oscillator for which 
good thermal averages are easily obtained and even with 
a separate transverse degree of freedom of each link, most 
of the work goes into sampling charge configurations. 

If this conclusion is correct, one would expect that the 
charge dynamics introduced in Section IIII Al in the pres- 
ence of thermal noise should be sufficient to recover elec- 
trostatic averages. Unfortunately the evolution of the 
electric field configuration is not a Markov process so 
further analysis is difficult. 



V. SUMMARY AND CONCLUSIONS 

We have made a study of a local Monte Carlo method 
developed by Rottler and Maggs for simulating sys- 
tems of charged particles. The unrefined method is lim- 
ited in the accuracy it can achieve, because the electric 
field is discretized onto a lattice. One must remove self- 
energy contributions due to charge spreading and correct 
the short-range interaction. All these can be done locally 
so the scaling of the method is O(N) where N is the num- 
ber of charges. 

1. Simulations of the 2D one-component plasma have 
been carried out in one of the first applications of the lo- 
cal off-lattice Monte Carlo method. With a short-range 
interaction correction extending a distance of about a 
tenth of the length L of each simulation cell from each 
charge, excellent results are obtained for the pair corre- 
lation function for T < 100. By making the method less 
local results can be obtained to lower temperatures. For 
example, with short-range correction a distance of L/5 
from each charge, good results can be obtained at least 
down to the transition temperature T ~ 140. 

2. Although there is a non- vanishing uniform electric 
field during our local simulation which is a consequence 
of the charge move update, we find that this produces 
no significant deviation from results obtained from an 
Ewald potential for which the dipole term is zero. It 
is not certain from this study whether this result will 
generalize to all systems containing many mobile charges 
but I suspect that the electric field will indeed average 
out because charge winding, which gives an electric field 
that cannot be inferred from the particle positions, is 
always important. 

3. In order to obtain accurate results at low tempera- 
tures (r ~ 140), a much larger rate of plaquette updates 
is necessary than at lower temperatures. This is corre- 
lated with the fluid-crystal transition. Plaquette updates 
are completed much more quickly than charge moves so 
this effect does not substantially alter the speed though 
of course charge autocorrelation times are larger at lower 
temperatures so simulations should be longer anyway. 
One might introduce a worm update to sample the trans- 
verse electric field more efficiently at low temperatures 
though it would seem to be unnecessary otherwise, at 



least for the 2D OCP. 

4. The charge autocorrelation times for the local 
method are much larger than for global methods. This 
is because there is, on average, a larger energy barrier to 
charge motion when Gauss' law is maintained locally. 

5. The autocorrelation times should be included in 
calculations of the effective speed. When this is done for 
system sizes that are typically simulated, we find that 
the gap between the local method and the 0(N 2 ) global 
method is considerably narrowed compared to the raw 
speeds, but that the local method remains significantly 
faster and the gap widens as the number of charges in- 
creases. We have not studied molecular dynamics meth- 
ods which in some cases have O(N) scaling. Prelimi- 
nary results by Pasichnyk and Diinweg [l6[ suggest that 
molecular dynamics adapted from the local Monte Carlo 
we have discussed here is a promising alternative to the 
P 3 M method. But, we expect Monte Carlo to have the 
advantage of relative ease of implementation and smaller 
scaling prefactor. The most significant contributions to 
the prefactor are charge assignment, local field updates 
and especially corrections to the lattice potential. Ex- 
cept at low temperatures, plaquette updates are not a 
significant factor in the speed. 

6. We also looked at the question of how the intro- 
duction of extra degrees of freedom can result in a faster 
method when naively we would expect the transverse de- 
grees of freedom should be allowed to relax across the 
simulation cell for each charge move. Charge mobility 
for intermediate temperatures is seen to increase sharply 
as the rate of plaquette updates is increased from zero. 
We have provided a simple model that, from a reason- 
able premise, reproduces this qualitative behavior. For 
high temperatures T < 2, the charge mobility is large 
even when independent plaquette moves are switched 
off. The fact that plaquette updates can be switched off 
at high temperatures strongly suggests that importance 
sampling of the simple transverse degrees of freedom can 
be accomplished entirely by charge moves, provided they 
are sufficiently mobile. 

7. We have introduced a continuum version of the 
charge move electric field update which allows us to see 
more rigorously that charge moves alter the transverse 
electric field. Following from the observation in the pre- 
vious point, an interesting problem is to carry out the 
average over electric fields using this dynamics alone and 
hence recover electrostatic averages. 

8. In an appendix, we have showed that a thermal 
average over Maxwell electrodynamics for a system of 
charges gives the same charge correlation functions one 
would obtain from electrostatics alone. As is well-known, 
the magnetic field has no effect on correlation functions. 
Also, the transverse part of the electric field averages out 
as in the original discussion by Maggs and Rossetto. 
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APPENDIX A: APPROXIMATION TO THE 
EWALD POTENTIAL IN TWO DIMENSIONS 



We now turn to the k space summation which, for small 
r is 

q 2 v—* exp(— 7r 2 k 2 /a 2 ) .„ , , 2 q 2 , , 9n 

The cross-terms in k • r vanish because the summation is 
over positive and negative integers. We can also remove 
a factor of r 2 because the terms in the summation are 
identical. Then we use the identity 



In this section, we find a good approximation to the 
Ewald potential in two dimensions for small charge sep- 
arations. Consider 2 charges q in a rectangular periodic 
cell with side lengths L x and L y . The whole system is 
charge neutral because a neutralizing charge is smeared 
uniformly over the cell. Leeuw and Perram [22| have 
found that the total energy can be expressed as 



a 2 A 



n i,j=l 
,2 „™/ ^2l,2 /„,2\ 



2ttA 



E 



exp(— 7r 2 k 2 /ar) 



k 2 



[1 — cos(27rk • 1-12)] 



-y(7 + ln^)-^ (Al) 

where the * means that we should not include the charge 
self energy in the summation over n = (n x ,n y ) and n, 
are integers. The area of the system is A = L x L y and 
7 is Euler's constant. The energy does not depend on 
constant a but the rate of convergence of the summa- 
tions does. This is the energy of the interaction between 
the charges and their periodic images, interaction of the 
charges with the background charge and the self-energy 
of the background. We would like to subtract off the 
background contributions. The interaction between two 
charges without images is just (q 2 / 2ir) ln(l / r) . So we 
compute 

q 2 

U bkg = l™U+—lnr 

because the interaction of the charges with their images 
are guaranteed to vanish in the limit. Hence the inter- 
action between the charges is given by Ui nt — U — £4fc g . 
We find 



2nA 



E 

k^O 



exp(— 7r 2 k 2 /a 2 ) 
k 2 



[-l + cos(27rk-r)] 



_1.( 7 + W) (A2) 

where the distance between the charges is denoted r. We 
approximate the exponential integral for small r 

£i(a 2 r 2 ) = -7 - 2 In or + a 2 r 2 + ... 



E2i 2 ! 2 X T 2 A 

exp— 7r k I a = 2 , ex P ~ a ^ n 



and approximate the exponential for large A, retaining 
the terms to order A. We find that the short range cor- 
rection energy Usrc, which approximates L^„ t , is 



2 99 
q* q'r z 



(A3) 



the terms depending on a cancel. The calculation in 
three dimensions is given in an appendix to 0]. 

We compare this result with the Ewald potential. The 
figure (|10p shows how the approximation to the short 
range correction given by equation (|A3[) differs from the 
Lekner interaction V from equation ([3]) as the charge 
separation increases. The measure we use is 

AV(x) - AUsrc(x) 



AV(x) 

where AUsrc{x) and AV(x) are energy differences as a 
charge steps outwards by half a lattice spacing from sep- 
aration x. The error is within 1% for charge separations 
of 1/5 of the cell width L x and for the highest accuracy 
low temperature simulations we performed, up to 2%. 



APPENDIX B: ELECTROSTATICS FROM 
THERMAL ELECTRODYNAMICS 

In this part we wish to demonstrate that the argument 
Maggs and Rottler have made for their rather unphysi- 
cal energy function holds also for full Maxwell electrody- 
namics. That is to say that even though the charges do 
not influence one another by instantaneous action at a 
distance, they nevertheless generate thermodynamic av- 
erages as though they do. This result has already been 
found by Alastuey and Appel [23| but I believe it would 
be useful to present the argument in this context. The 
main difficulty is in finding an energy functional that 
includes the constraints that depend on the static con- 
figuration of charged particles. 

The charges are not subject to an external field but are 
influenced by other charges in our hypothetical box in a 
heat bath. Since the fields are dynamical quantities we 
must include them in our discussion. So, the Lagrangian 
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FIG. 10: Plot of the normalized error between energy dif- 
ferences calculated using the Lekner potential and the short 
range interaction potential (|A3[) . We have used parame- 
ters appropriate to the simulations we have performed, with 
32 x 32 square plaquettes covering the cell and lattice spacing 
a — 0.342. The energy differences are for charge steps in the x 
direction through a distance of one-half of the lattice spacing. 
L x is the width of the periodic cell in the x direction. 



for a system of N relativistic charges qi for i = 1,2, ... ,N 
interacting electrodynamically is 



L = 



N 

£ 

i=l 



-m , 



1 - v? + J d 3 r (-If^F^ + J„A^ 



where {v^} are the particle velocities and {men} their 
rest masses. F^ = d^A v — dvA^ is the electromagnetic 
field tensor defined in terms of 4- vector potential A^. 
The 4-vector current has components J° which is the 
charge density and spatial components in 3-vector J = 
J2i <K r — Yi)qiVi. We have taken c = 1. 

The form of the Lagrangian ensures that the homoge- 
neous Maxwell equations are obeyed [2(| 

{ v x E = — at 

The first step in getting a Hamiltonian is to find the con- 
jugate momenta from the Lagrangian. For the charges, 
the momenta are 

dL m 0i w l 
Pi = ^— = , 9 + %Aj 
ovi ^/l - v f 

where = A(rj) and for the fields 



IP 



dC 
doA, 



-F, 



Ofi 



where C is the Lagrangian density. We notice that II 
vanishes because F^ u is antisymmetric. The equations 
of motion are 

d u F^ + J v = 0. 



If we go straight ahead and find the Hamiltonian for 
the fields and the interaction term 



H F = I d 3 rn • A - I d d r[ —F^F^ + J„A» 



we get 



H F = I d 3 r ( ^ (B 2 + E 2 ) - J • A 



in terms of the electric and magnetic fields having used 
an equation of motion after partial integration to write 

n A = FoiFoi + A J°. 

In this form, we cannot perform an average over the 
electric field because it is constrained by Gauss' law 

V • E = J° 

which is formally an equation of motion although it is 
a static constraint on the field. We shall now render 
the Hamiltonian in a form that brings out the important 
constraints. To do this, we first notice that the electric 
field can be split into rotational and irrotational parts, 
exactly as we did in the introduction. Using the same 
notation as before 



E — En 



E, 



We also recall that the Lagrangian for electrodynamics 
is gauge invariant. We can fix the gauge to reduce the 
arbitrariness in A. For our purposes it is useful to choose 
the so-called Coulomb gauge V • A = 0. This can always 
be done. The reason for choosing this gauge is that the 
splitting of the electric field into transverse and longi- 
tudinal parts is mirrored in the vector potential in the 
following nice way 



En = -VA° Ej_ — — 



dA 
~dt' 



We see that Gauss' law in terms of Aq is Poisson's equa- 
tion which has the solution 



A fat) = J d 3 r 



3 , J°(r',t) 



4tt r - r' 



which is the instantaneous Coulomb interaction. The 
transverse field is independent of the Gauss' law con- 
straint. It is also significant that, as before, the electric 
field energy splits into independent transverse and longi- 
tudinal energies 



d 3 r E 



d d r El 



Udw J °^V), 

2 / 47r r-r' 



Henceforth we denote the longitudinal field energy by V . 
Now that we have put the longitudinal dependence into 
A we would like to have a momentum variable that is in- 
dependent of A . We choose = II — WA° which must 



17 



satisfy the constraint V • 11^ = that follows directly 
from the gauge constraint. We must get the Hamilto- 
nian from this new momentum. For the field quantities 



d v IT • A - L 



Hp 



where Lp is the Lagrangian excluding the field- 
independent parts, so that the complete energy function 
(including the free particle and interaction terms) is 



d 3 r 



Ini + i(vx Af + v 



The partition function in the canonical ensemble with 
inverse temperature /3 is 

i i 

x / PA(5(V • A) exp(-fiH) 



and the limits of the particle momentum integrations 
are infinite. The field momentum integration factorizes 
from the rest. Then we make the change of variables 
Pi = P« — 9iAj for each particle. The Jacobian for this 
transformation is unity because the vector potential docs 
not depend on the momenta. So, the integral over the 
vector potential also factorizes and we are left with 



Z(P) = Z(P) P Z(P) n Z({3) A / Hd 3 ^ cxp(-pV) 



which, up to a configuration independent factor is the 
electrostatic partition function. This is essentially a gen- 
eralized version of the well-known Bohr- van Leeuwen the- 
orem (see, for example [HI)- 
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